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Abstract 

This paper develops efficient ensemble Kalman filter (EnKF) im¬ 
plementations based on shrinkage covariance estimation. The forecast 
ensemble members at each step are used to estimate the background 
error covariance matrix via the Rao-Blackwell Ledoit and Wolf esti¬ 
mator, which has been developed specifically developed to approxi¬ 
mate high-dimensional covariance matrices using a small number of 
samples. Additional samples are taken from the normal distribution 
described by the background ensemble mean and the estimated back¬ 
ground covariance matrix in order to increase the size of the ensemble 
and reduce the sampling error of the filter. This increase in the size 
of the ensemble is obtained without running the forward model. Af¬ 
ter the assimilation step, the additional samples are discarded and 
only the initial members are propagated. Two implementations are 
considered. In the EnKF Full-Space (EnKF-FS) approach the as¬ 
similation process is performed in the model space, while the EnKF 
Reduce-Space (EnKF-RS) formulation performs the analysis in the 
subspace spanned by the ensemble members. Numerical experiments 
carried out with a quasi-geostrophic model show that the proposed 
implementations outperform current methods such as the traditional 
EnKF formulation, square root filters, and inflation-free EnKF imple¬ 
mentations. The proposed implementations provide good results with 
small ensemble sizes (~ 10) and small percentages of observed com¬ 
ponents from the vector state. These results are similar (and in some 
cases better) to traditional methods using large ensemble sizes (~ 80) 
and large percentages of observed components. The computational 
times of the new implementations remain reasonably low. 

Keywords: EnKF, shrinkage covariance estimation, backgronnd errors, 
sqnare root filter 

1 Introduction 

Seqnential data assimilation estimates the cnrrent nnknown state G 

of a physical system as follows. The background (prior) state G R”^^ 
is given by a physical model initialized with the best estimate of the state at 
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a previous time: 


X — Xcurrent — -^tprevious^icurrent (Xprevious) • (1) 

The background errors are assumed to be unbiased and normally distributed 

^ = x'-x'™'^eAr(0„, B) (2) 

where n is the dimension of the model state, 0^ is the vector in the q- 
dimensional space whose components are all zeros, and B G is the 

unknown background error covariance matrix. Observations of the true state 
are taken 


y = H + e e , (3) 

where Ti : IR"'^^ —)■ IR™^^ is the observation operator, and m is the number 
of observed components. The observational errors are assumed to be nor¬ 
mally distribnted e ~ M {Om, R-), where R G is the observation error 

covariance matrix that is assnmed to possess a simple strnctnre (e.g., block 
diagonal) and therefore its inverse can be easily compnted. 

Under the Gaussian assnmption on data and model errors the negative 
logarithms of the a posteriori probability density is the 3D-Var cost fnnction 

EH: 

>7 W = ) ■ ||x - x*|[ 3 _, + t ■ ||y - W (x) IIh -1 . (4) 

The maximum likelihood estimate of the state is obtained by minimizing the 
cost fnnction Q, i.e., the analysis state x“ G IR"^^ is the solntion of the 
following optimization problem: 

x“ = argmin J (x) . (5) 

X 

The solntion of (|^ over the snbspace spanned by the ensemble members is: 

X“ = x"- + , (6) 

where = B^/^ is a set of basis vectors satisfying = B, and the vector 

of weights is given by 

= Vg ■ (R + Vb ■ Vg) ■ (y - , (7) 
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where Vb = H - G and Yl = 1-1'. By propagating in time the analysis 

(|^ a new prior (background) state is obtained for the next assimilation cycle. 

Since B is unknown the direct use of (§ is infeasible. Methods for esti¬ 
mating the background error covariance matrix have been developed [T3l |2l] . 
However, in the absence of prior information about the true structure of B 
biased estimators are often obtained. Buhener [7] discusses the impact of 
having biased estimators for the background errors in the assimilation pro¬ 
cess. Ain the context of sequential data assimilation an ensemble of model 
realizations [18] is built in order to represent the background error statistics: 

x'' = [x;, x‘, ...,xy (8) 

where G is the i-th ensemble member (model run) in the background 
stage. Estimates of and B are obtained via the empirical moments of the 
ensemble ([^ 


N 




nx 1 


(9) 


Z=1 


and 


B ^ = S • G R" 


( 10 ) 


where the matrix of scaled member deviations is 

1 


S = 




[X' 


6 


1?;] e R' 


nxN 


( 11 ) 


where Ig is the vector in the g-th dimensional space whose components are 
all ones. The set of basis vectors (11) does not span the full space of model 
errors. By replacing the estimators and (10) in (§, the analysis can be 
approximated as follows: 


x“ = xF S ■ a G R' 


nx 1 


( 12 ) 


where 


CK = V ■ (R + V ■ V^) ^ ■ (y - (x^)) G R 


)Wxl 


(13) 


and V = H ■ S G The analysis state (12) is the optimal solution 

in the space spanned by the ensemble members but in general it is not 
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optimal in the model space. Since n 3> iV information contained in Q is 
not represented by the set of basis vectors S. From another point of view, 
the ensemble covariance matrix is rank dehcient and therefore there are 
insufficient degrees of freedom to explain the full error. This problem can 
be alleviate by making use of localization techniques ilElES]. However, 
since the structure of B remains unknown, the use of localization generally 
increases the bias in the background error estimation. 

There is an opportunity to avoid the intrinsic need of inflation in ensemble 
based methods by replacing the covariance matrix (10) with a more accurate 
and well-conditioned estimate of B. We do not want to impose any kind of 
structure on B since it will make our approach sensitive to problems faced by 
current implementations. Instead, we seek to capture most of the informa¬ 
tion contained in the ensemble. Shrinkage covariance estimators developed 
to estimate high-dimensional covariance matrices from a small number of 
samples |30] £t very well in the context of sequential data assimilation. 

The remaining part of the paper is organized as follows. Section re¬ 
views ensemble based data assimilation and shrinkage covariance estimation. 
In section the two novel implementations of the ensemble Kalman hlter 
based on shrinkage covariance estimation are proposed. Experimental re¬ 
sults making use of a quasi-geostrophic model are given in section Section 
|5] summarizes the conclusions of this work. 


2 Background 

In this section we review relevant concepts with regard to shrinkage covari¬ 
ance estimation and ensemble based methods in sequential data assimilation. 

2.1 Covariance estimation 

Many problems in science and engineering require an estimate of a covariance 
matrix and/or its inverse, where the matrix dimension n is large compared to 
the sample size N. Different applications ranging from variational [HEl] to 
sequential data assimilation rely on accurately estimated covariance 

matrices. 

Let {si, S 2 , ..., sw} be a sample of independent identical distributed n- 
dimensional Gaussian vectors 

Q) 
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A common approach is to estimate Q G by the sample covariance 

matrix Cg 


C. = 


N-1 


N 

i=l 


si e 


(14) 


Cs is the maximum likelihood estimator when it is invertible (SO]. However, 
under the condition n S> A^, this is not the case. The simpler thing to do 
in order to deal with the rank-dehciency of C* is to impose some structure 
(i.e., localization in ensemble based methods). However, in the absence of 
prior information about the true structure of Q, C* will poorly describe the 
correlations between different components of the samples {sj}i<j<Ar. In or¬ 
der to improve estimation of covariance matrices many methods have been 
proposed in the literature based on tapering procedures | 8 ], [ 10 ], minimiz¬ 
ing the log-determinant divergence |36], and greedy methods ^7\. Another 
class of well-conditioned estimators is based on shrinkage approximations 
la [151 USmUEDlES]. These approximations express the estimated covari¬ 
ance matrix as a weighted average of some target matrix T G R"'^"' and the 


empirical covariance matrix (14). To better understand this assume that the 


components of Si are uncorrelated for 1 < i < iV. A simple estimate of Q is 
given by 


T = 


tr (Cs) 


I 


n 


nxn 5 


where I„xn is the identity matrix in the n-dimensional space. Note that 
this structure will reduce the variance but will increase the bias when the 
diagonal assumption is not fulhlled. A reasonable trade-off is achieved by 
the shrinkage of Cs towards T and provides the followng class of estimators 


C = 7 ■ T (1 - 7) ■ C, G R" 


(15) 


where 7 G [0, 1]. The problem is then reduced to hnd an optimal value for 7 
in which the squared loss 


E 


C-Q 


(16) 


is minimized, where ll^llp denotes the Frobenius norm. There are many 
shrinkage based estimators derived from the minimization of (16) subject to 
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( p!5| . We restrict our exploration to three well-accepted methods: the Ledoit 
and Wolf estimator EDI. the Rao-Blackwell Ledoit and Wolf estimator m 
and the oracle approximating shrinkage estimator mn 

The distribution-free Ledoit and Wolf (LW) estimator [30] has been proven 
more accurate than the sample covariance matrix and some estimators pro¬ 
posed in hnite sample decision theory. Moreover, it is better conditioned 
than the true covariance matrix PI. The optimal 7 value proposed by this 
estimator is 


jLiv = min 


Eii ||c. 


T\ 

si 


iV2 . 





(17) 


and the LW estimator C^w is obtained by using 'Jlw in 

The Rao-Blackwell Ledoit and Wolf (RBLW) estimator [TTl [T2] provably 
improves the LW method under Gaussian assumptions. The motivation of 
this estimator is that, under Gaussian assumptions, all the information re¬ 
quired in order to get a well-conditioned estimate of Q is contained in C^. 
The proposed value for 7 is 


IRELW = mm 


N-2 

n 


■tr(C2)+tr2(a) 


(AT+ 2)- tr(C2)-h:^ 


(18) 


and the corresponding estimator Crblw is obtained by replacing (18) in the 
equation (15). In addition, in [T^ Theorem 2], it is proven that 




2 ■ 



2 ■ 

E 

CrbLW ~ Q 

F 

< E 

Clw ~ Q 

F 


which rigorously shows the RBLW estimator to be a better approximation 
of Q than the LW estimator under the Gaussian assumption. 

The oracle approximating shrinkage (GAS) estimator |TT] is an iterative 
approximation of the unimplementable oracle method m Section 3]. The 
optimal 7 at each iteration j is given by 


Cj — 7 j • T -t- (1 7 j) ■ Cs , 

(1 - i) ■ f {C, ■ C.) + tr^ (7) 


7i + l 


(AT + 1 - I) . tr (C, . C.) + (1 - f) ■ te (C,) 


(19a) 

(19b) 
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where the initial estimator Cq can be any estimator (i.e., the LW, the RBLW, 
or even the sample covariance matrix). 

All the estimators presented in this section provide well-conditioned ap¬ 
proximations to the unknown covariance matrix Q. We center our attention 
on the RBLW estimator since in high dimensional problems, such those found 
in data assimilation, this estimator can be implemented easily, and under the 
Gaussian assumption it provides better approximations than the LW estima¬ 
tor. 


2.2 Sequential data assimilation methods 


Ensemble data assimilation methods are widely used in applications to weather, 
oceanography, and climatology 021 . These methods represent the back¬ 
ground error statistics by the empirical moments of the ensemble (|^ and 
(10)[22]- The trajectory of each ensemble member and the dispersion of the 


overall ensemble around the background state provide meaningful informa¬ 
tion about the background error distribution. One of the most important 
advantages of ensemble DA is the flow-dependent forecast error covariance 


matrix [7]. When observations are available the analysis state (12) is com¬ 
puted. The next step is to generate an ensemble which describes the analysis 
uncertainty around this optimal state. We briefly discuss three implementa¬ 
tions that achieve this, the ensemble square root hlter, the ensemble trans¬ 
form Kalman hlter, and the (basic) ensemble Kalman hlter. 

In ensemble square root hlters [39], the new analysis ensemble in the 
analysis is built as follows: 


X“ = x“ ® 1^ + S 


NxN 


■ (R + V ■ V^) ^ ■ V 


1/2 


G 


nxN 


( 20 ) 


As expected the analysis ensemble members live in the subspace space spanned 
by the background ensemble, this is, the space spanned by the columns of 
S. All possible information that can be obtained from the model states is 
contained in this set of basis vectors. 

The covariance matrix in the observation space is 

Wobs = R + V ■ G R™^”^, (21) 

the linear system Wobs ■ Zv = V G R™^^ can be solved via the iterative 






Sherman Morrison formula (ISMF) [32] 


hW = fl + vj.ui'-'' 


-1 


u 


(k-1) 


e 


z(fc) = z(fc-i) _ hW . [v^ . G E 

U(fc) = ijik-i) _ hW . ^ 


mxN 

mxN 


(22a) 

(22b) 

(22c) 


for 1 < A; < iV, where Z^o) = ■ (y - H ■ x^) and = R"^ ■ V. The 

implementation ( [22| ) requires no more than O {N'^ ■ m) long computations, 
where and are the fc-th column of the matrices V and 

respectively. Moreover, by applying the singular value decomposition to 


• Zv = Uz ■ Sz ■ U4 G E 


NxN 


G E 


NxN 


(23) 


the square root in (20) 


r = Uz-[l7VxiV-Sz]'/'-Ui, 

can be efficiently computed with no more than O {N^) long computations. 

Another widely used square root hlter implementation is the ensemble 
transform Kalman hlter (EnTKF) [33|. By making use of the matrix identity 


I - V ■ [R + V ■ V'^] ^ = [I + V • R ■ V"^] ^ G E*^ 

and the singular value decomposition 


= Uv • Sv ■ G E^^’” 


the analysis state (12) can be written as follows 

x“ = x^ + S ■ /3 

with the optimal weights (3 G E^^^ given by 

/3 = UvSv (livxiv + Sv-Sv)''-V^-yR. (y - H ■ x') . 
The new ensemble members are built as follows 


X“ = x“ 0 + S ■ Uv ■ 


NxN 


+ Sv ■ Sy- 


1/2 


■U^ 


V- 


(24) 


(25) 


(26) 


(27) 


(28) 
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In the ensemble Kalman filter (EnKF) [T7j each ensemble member (|^ and 
the observations (|^ are treated as random variables, and the i-th. ensemble 
member is npdated as follows: 

x“ = x‘ + S-/3,eR”'<‘, » = (29) 

where the optimal weights for i-th ensemble member are given by 

/3, = ■ (R + V ■ V^) - H ■ G (30) 


and 


y‘~M (y, R) 


(31) 


The addition of the pertnrbed observations (31) in the analysis provides 


asymptotically correct analysis-error covariance estimates for large ensem¬ 
ble sizes and makes the formulation of the EnKF statistical consistent Ha. 
However, it also has been proven that the inclusion of perturbed observations 
introduces sampling errors in the assimilation [21121] • 

One of the important problems faced by current ensemble based methods 
is hlter divergence due to the insufficient degrees of freedom (A^ n). To 

alleviate this dehciency localization is used to impose structure on the sample 


covariance matrix (10) according to the physics of the model. Intuitively, the 


correlations between individual model variable errors decays with distance, 
e.g., exponentially: 


Pij = exp 


2 ■ L2 


(32) 


where d{i,j) is the physical distance between the locations of the i-th and 
j-th model components, and L is the localization radius. The localized co- 
variance matrix is obtained as 

= poP^ e (33) 

where o is the Schur product. This method can be impractical since it relies 
on the explicit computation of P^. Moreover, there are no guaranties that 
this method captures the true structure of B. More sophisticated covari¬ 
ance estimation methods have been proposed in the context of data assim¬ 
ilation. A classic approximation is the Hollingworth and Lonnberg method 
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|2S] in which the difference between observations and background states are 
treated as a combination of background and observations errors. However, 
this method provides statistics of background errors in observation spaces, 
requires uniform observing network (not the case in practice), and the re¬ 
sulting statistics are biased towards data-dense areas. Another method has 
been proposed by Benedetti and Fisher [3] based on forecast differences in 
which the spatial correlations of background errors are assumed to be similar 
at 24 and 48 hours forecasts. This method can be efficiently implemented in 
practice, however, it does not perform well in data-sparse regions, and the 
statistics provided are a mixture of analysis and background errors. Since 
the structure of B remains unknown, assumptions made about its structure 
may increase the bias in the estimate. Furthermore, the balance of variables 
in the model with some physical meaning can be disturbed when localization 
is utilized [34] . 

A different approach is based on the 3D-Var cost function in the ensemble 
space. Any vector x G in the ensemble subspace can be written as 

X = + U ■ CK (34) 

where U is the matrix of anomalies 


U = [x5 — X^, X2 — x^, ..., x^ — x^] G R 


nxN 


(35) 


and a G R'^^^ is a vector to be determined. The columns of U and S span 
the same space. Using (34) the 3D-Var cost function (|^ in the ensemble 
space reads 


Jens («=) = ^ ■ l|U ■ Q:||b-i + \ ■ IItQ ' «IIr-i , 
where the optimal value of the control variable cx. 

a' ^ aTgmmJgr^(a) , 


(36) 

(37) 


provides the analysis state in (34) 


x“ = x^ 


U a* 


(38) 


Two recent formulations based on this approximation are the finite size 
anddual ensemble Kalman filters [6]. These formulation avoid the intrin¬ 
sic needed of inflation by choosing Jeffrey’s prior for background errors: 

iP(x',B) =iPj(x') ■iPj(B) , 
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where the parameters and B are assumed to be independent. 

In the case of the hnite size ensemble Kalman hlter (EnKF-N) the cost 
function in the ensemble space reads 

(«) = ^ lly - + U • a) ||^_1 + y ■ log (^1 + ^ + llaf ^ . (39) 

Minimization of this cost function provides the optimal weights in the en¬ 
semble space 


a" = arg min J™, (a) , 


(40) 


and the analysis is computed via (38). The projection of the analysis co- 
variance matrix on the ensemble space is approximated by the inverse of the 

(41) 


Hessian of (39) at the optimal value (40). The Hessian reads: 

iT 

■ jx ■ r 

|2 


= [H-Ur-R-'-H-U 


_l_ _ (l + jv + I|Q^II ) • lyfxAf 2 ■ OL ■ OL^ ^ 


(l + yr + I|Q:||^) 
The analysis ensemble is generated as follows: 


X“ = x“ ® l)v + U ■ 


_n 1/2 


-1)- I (42) 


where $ G is an arbitrary orthogonal matrix which preserves the 

ensemble mean ([^. 

Another approach is based on the dual formulation of the cost function 

® [ 6 ] 

I>=” (C) = l|y - H . x‘|P + c • (l + 4) + A' ■ log (?) - A'. (43) 


IW^-^ ’ ^ y ' Ny 
where the weighted covariance matrix reads 

= R + i • V ■ G R™' 


C 


C 


The dual optimization problem is one-dimansional 

C = arg min subject to C ^ ( 0) 


N 


1 + iv 


(44) 


(45) 
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The optimal state is computed as follows: 

x“ = x"' + U • ■ R-' ■ V + C* ■ 1nxn]~' ■ • R-^ ■ [y - H ■ x^] , (46) 

and the following analysis ensemble is built: 

X“ = x“ 0 1^ + U ■ |(iV - 1) • [V^ ■ R-i ■ V + C ■ (47) 

In this paper we consider a different representation of the background 
error statistics by making use of shrinkage covariance estimation. The idea 
is not to impose any structure on but to obtain a well-conditioned es¬ 
timator B of the background error covariance matrix B wherein using all 
the possible information brought from the ensemble members. Samples from 
the distribution Af ^x^, B^ are taken in order to better represent the error 
statistics and to increase the number of degrees of freedom. Two novel EnKF 
implementations based on the Rao-Blackwell Ledoit and Wolf estimator are 
presented in the next section. 


3 Ensemble Filters Based on Shrinkage Co- 
variance Estimators 

In this section, we propose two efficient implementations of the EnKF based 


on the RBLW estimator (18). As mentioned before, we do not impose any 
kind of structure on since the information brought by the ensemble mem¬ 
bers is more than only background errors 

P'’ = B + Q + Cg , 

where Q is the covariance of model errors and C G R"^"" is the covariance 
matrix of additional errors whose sources are unknown for us. Errors coming 
from different sources are assumed to be uncorrelated. We seek to exploit the 
information brought by ensemble members and use the RBLW covariance es¬ 


timator (18) to build a covariance matrix that captures all error correlations. 


The standard form of this estimatordepends on the explicit representation of 
P^. The efficient implementation for high-dimensional covariance matrices 


presented in section 3.1 avoids the explicit computation of P^. Section 3.2 
discusses two EnKF implementations based on the RBLW estimator. Sec¬ 
tion [3^ develops an efficient sampling method in high dimensions for drawing 


13 








samples from the prior error distribution based on the RBLW estimate. Fi¬ 


nally, section |3.4| discusses the similarities and differences between the two 
proposed implementations. 

3.1 RBLW estimator for covariance matrices in high- 
dimensions 


Consider the sample covariance matrix (|10|). In the context of data assimi- 

, (48a) 


lation the RBLW estimator (15),(18) reads 

B = 7b ■ • I + (1 — 7b) ■ ^ K- 

where 

tr (P'’) 


/^B ~ 


n 


( 


7b = mm 


iV-2 


tr +tr2 (P^) 


(AT+ 2 ) 


tr ( [P^l^ 


tr^ 


(p*>) 




/ 


(48b) 

(48c) 


Since the dimension of the model state is high (n ~ 0(10^)), the direct 


computation of (52) is impractical as it requires the explicit representation 
of the sample covariance matrix P^. An alternative manner to compute 
tr (P^) and tr ([P^]^) is proposed. Consider the eigenvalue decomposition of 


P^ = Up. ■ Sp. ■ G : 


(49) 


where Sp. G R"^" is a diagonal matrix whose diagonal components cTj, for 
1 < i < n, are the eigenvalues of P^ and Up. G R”^"" is a set of orthogonal 
basis vectors spanning the ensemble space (since P^ is rank deficient). By 

dehnitiont tr (P^) = Since there are only 

N — 1 eigenvalues different from zero we obtain: 


N-l 


tr 


2 = 1 


(p") = ‘d[p‘l) = 


Af-l 


i=l 
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and the computations in the set of equations (48) can be efficiently performed 
whenever the hrst iV — 1 eigenvalues of can be easily obtained. Consider 


the singular value decomposition (SVD) for the set of basis vectors 01 


S = Us ■ Ss ■ Vi e R 


nx N 


(50) 


where Ss G is a diagonal matrix holding the singular values a* of S, 

for 1 < i Likewise, Us G and Vs G are the left and right 

singular vectors, respectively. Since = S ■ we have Spt = Ss ■ Sg and 

N-l N-1 


tr(P‘) = 


CTi 


i=l 

N-l 


i=l 

N-l 


tr ( [P'']' 




(^i 


2=1 


2=1 


The computational effort of the SVD decomposition (50) is (9 {N^ ■ n). The 


traces in (48) can be computing without calculating the sample covariance 
matrix P^ by making use of the inexpensive SVD decomposition of S. None 
of the singular vector of S are required, but only the singular values a*, for 

(51a) 

(51b) 


1 < i < N — 1. The parameter values in (48) are computed as follows: 

•sr-^N—l -^2 




n 


( 


7b = mm 


V' 


N-2 s^N-l ;-4 1 

sr-^N—l -^2 

2 

(iV + 2) ■ 

1-1 

M 

■^N-1 -2 

-12 " 

n 



, 1 


/ 


With = /ig ■ 7g and h = 1 — yg the estimated covariance matrix 
B = (^ ■ Ux„ + 5 ■ S ■ G . 


IS 

(52) 


3.2 EnKF implementations based on the RBLW esti¬ 
mator 


By replacing the estimated error covariance matrix (52) in (29), the EnKF 
analysis in matrix form becomes 

-1 


X“ = X'^ + B R + H B H 


D G E 


nxN 


(63) 
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where the matrix of innovations D G R’ 


mxN 


IS 




( 54 ) 


and the data y* for 1 < i < is drawn from the distribution (31). We have 

■ (R + H-(yp-I„xn + 5-S-S^) 

X“ = X^’ + E-n-Zg + yp-H^-Zg, (55) 


where E = ^5 ■ S G H = H ■ E G and Zg G R^ 

the solution of the linear system 


(r + n-n^)-ZB = d, 

r = R + (y9-H-H^ G R”^^™ 


mxAT is giygB by 

(56) 


When H possesses a simple structure (e.g., indexes to observed components 
from vector states) the matrix E also has a simple structure (since in practice 
R is block diagonal). By letting • El G R™^^ and Zg^ = E^^ ■ D, 


the linear system (56) can be efficiently solved via the ISMF with no more 


than O ■ m) long computations. When E has no special structure its 
inverse can be calculated off-line. 

In order to obtain a better representation of the background error statis¬ 
tics (uncertainty) about the background state (|^ additional samples can be 
taken from the distribution 




N (x^ B) G R^ 


) nx 1 


1 < i < K. 


(57) 


This yields to a new ensemble formed of two kinds of members, real and 
synthetic. The real members are obtained by model propagation of 

the previous analysis ensemble. The synthetic members {x('}^^ are artifi¬ 


cially built by taking samples from the distribution (57) and do not require 
additional model runs. 

The artificial increase in the size of the ensemble is therefore a relatively 
inexpensive modality to bring in additional degrees of freedom in the so¬ 
lution of the optimization problem ([^. Figure exemplifies the effect of 
additional members using two-dimensional projections of ensemble member 


states from the Lorenz 96 model. Figure la shows the spread of the real 
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ensemble members (for the background uncertainty around x°). Figures lb 


shows the distribution when artihcial members are added to the background 
ensemble, resulting in a better representation of the background error and 
therefore a decrease in the sampling error. 


“T-1-1-1-:- 

¥ 

* 

* * * 

**1*^^** ** ^ 

* « ♦ ** ** 

* * 

* * % 

¥ 


(a) K = 0 


++ 


-+ ^ 

5*++^ 


+ 




+ 


* + +%+t+* 

+#* + % ++ + 
* + 


(b) K = 120 


Figure 1: Error distribution for different values of K using two-dimensional 
projections of ensemble members from the Lorenz 96 model. The number of 
real members is iV = 40. In the plots * represents real ensemble members, 
-|- artihcial members, and o is the ensemble mean. 


The analysis state is now computed in the subspace spanned by both real 
and artihcial members: 


x“ G span {x?, x^, ..., x^, x?, x^, ..., x^} . 
The extended ensemble reads: 

^ — [X;^, X2, ..., x^, X;^, X2, ..., x^j t ja 


(58) 


(59) 


where Nk = N+K. Similar to (55), the background ensemble (59) is updated 
using 


X“ = X'' + E ■ n ■ Zg + (p ■ • Zb G IR' 


nxN 


(60) 
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where 


E = 

■ S G , 

(61a) 

n = 

H ■ E G , 

(61b) 

Zg = 

(^r + n ■ ^ d g r™^^ , 

(61c) 

s = 

X - x^ 0 G R”^^'= . 

(61d) 


The solution of (61c) can be efficiently obtained via the ISMF. Even if the 


artihcial members (57) are used in the covariance approximation for the 


analysis step, according to (60) only the states of real members are adjusted. 


After the assimilation step the artihcial members are discarded and only the 
real members form the analysis ensemble. This strategy does not increase the 
number of ensemble members to be propagated by model runs (and there¬ 
fore, the computational effort). When more computational resources become 
available, or when some real members are lost due to hardware failures, se¬ 
lected artihcial members can be updated and propagated as well. In this 
case they become real members during the next assimilation step. Moreover, 
some real members can be replaced by artihcial members such as to refresh 
the ensemble directions, e.g., in order to prevent hlter divergence. 

To summarize the above discussion, the implementation of the ensemble 
Kalman hlter based on the RBWL estimator with weighted covariance matrix 
in the observation space (EnKF-FS) consists of the following steps: 


Estimate the background error covariance matrix (52) based on the 
samples for 1 < i < A^. 


2. Compute the innovation matrix D according to (54). 


3. Draw K artihcial members from the distribution (57). 


4. Compute the set of matrices (61). 


5. Perform the assimilation (60). 


6 . Propagate the ensemble members 


xLxt = M 


^current—^^next V'^CUrrent/ ’ 


until the next assimilation step, for 1 < i < iV. 
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Another efficient implementation of the hlter based on the RBLW esti¬ 
mator can be obtained via the 3D-Var cost function in the ensemble space 


(36). This involves projecting the weighted covariance matrix (21) onto the 


ensemble space (58) rather than onto the observation space. 


After taking the samples (57) and building the ensemble (59) a new set 


of basis vectors can be built as follows 


U= [x? 




~b _ -b 

X^ X , Xi X 


5 ^1 




- X^] e , 


(62) 


where x^ is given in ([^. By replacing (62) and the RBLW estimator of the 
background error covariance matrix (52) in the 3D-Var cost function (36) we 
obtain 


jTens (A) — 2 


U A 


2 1 
B-I 2 ’ 


D-Q A 


R-i 


(63) 


where A G is the matrix of weights whose i-th column represent the 


coordinates of the i-th ensemble member in the space (58), for 1 < i < N^, 


and Q = H ■ U G The resulting 3D-Var optimization problem is 


A* = arg min Jens (A) , 

A 


(64) 


and has the solution 


A* = 


■ Zgu + R-^ ■ Q R ^ D G R^'= 


-1 


.xN 


(65) 


where Zg^ = B ^ ■ U G The resulting analysis ensemble is 


X“ = X'’ + U ■ A* G R' 


nxN 


( 66 ) 


To summarize, the implementation of the ensemble Kalman hlter based 
on covariance estimation with weighted covariance matrix in the ensemble 
space (EnKF-FS) consists of the following steps: 


1. Estimate the background error covariance matrix (52) based on the 
samples |x^W}^_^, for 1 < i < A^. 


2. Compute the innovation matrix D according to (54). 


3. Draw K artihcial members according to (57). 
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4. Compute the matrix of optimal weights (65). 


5. Perform the assimilation (66). 


6 . Propagate the ensemble members 


= Mt 


until the next assimilation step, for 1 < i < iV. 


3.3 Sampling in high-dimensions based on the RBLW 
estimator 


Both implementations discussed in Section |3.2 use samples from the distri¬ 
bution (57). Such samples can be generated as follows: 


SJ = x‘ + B'/" = I.„„ + i . s . S’’)''" ■ 


(67) 


for 1 < i < K, where ~ A/” (0„, Inxn)- However, this computation requires 
the explicit representation in memory of the estimated error covariance ma¬ 
trix B, which is prohibitive for high-resolution models. Moreover, the square 
root matrix is required making the use of (67) impractical. 


We need an equivalent strategy to obtain the samples (57) that requires 


a reasonable computational effort and does not use a full representation of 
the covariance matrix B. Toward this end consider the random vectors 


^2 


A7(0„, 

A^(0]v, Inxn) G 


and let 


COV (g>^f = OnxN , 

COV (^2, ^l) = 8) ^ = Ojvxn • 


We make the following substitution in 

B''" ■ ~ vV ■ «.'+ v'S ■ S ■ «?. 
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This does not change the statistics since 


= v?- +vV^- 

Cov(4l,4.')=IuXn Cov(€l42)=0„xiV Cov(€2,4l)=Ojv 

+<5 ■ S ■ ■S'^ = (^ • I„xn + 5 ■ s ■ S'^ 


= B. 


The artificial ensemble members are obtained as follows: 

X^ = x' + y^-I„xn-^- +V^-S-^,^ t = l,...,K. (68) 

The components of the random variables and are drawn indepen¬ 
dently from the standard normal distribution A/" (0, 1). For large model reso¬ 
lutions the components of the random vectors can be prepared independently 
taking advantage of parallel computations. Moreover, the random vectors 
and can be sampled prior the assimilation process in an off-line computa¬ 
tion. 

The estimated error covariance matrix is never represented explicitly in 
memory. Instead, the estimator B is represented via the triplet 

B = [ip, /i, S] . 


which contains two scalars and one matrix of dimension n x N. In addition, 
the scalars ip and p are computed making use only of the matrix S. This 
data is sufficient for correct sampling from the distribution (57). 


3.4 Comparison of EnKF-FS and EnKF-RS versions 
of the filter 

Although both EnKF-FS and EnKF-RS methods are based on the EnKF 
equations and RBLW estimator, their underlying theoretical properties are 
slightly different. To facilitate the comparison of the two proposed imple¬ 


mentations we bring the EnKF-FS analysis equation (60) to the form 


x“ = + b^/ 2 . cKg e 


nx N 


(69) 
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where the weights CKg G are given by 


CKg = 


R + H B H' 


-1 


D G 


nxN 


It is readily apparent from equations (69) and (66) that EnKF-FS and EnKF- 


RS implementations differ in the number of degrees of freedom used in the 
assimilation process. In the EnKF-FS approach the columns of B^/^ serve 
as the basis set for generating an ensemble of background deviations. Since 
the estimated background error covariance matrix is full-rank the optimal 


solution (46) is searched for in the full space. The matrix identity 


B 1 R-^ H 


1 -1 


R ^ = B H 


R + H B 


-1 


together with (69) reveal that the weighted covariance matrix of the EnKF- 


FS implementation 


W = B"^ + R 


H G E*^ 


(70) 


is related to the weighted covariance matrix of the EnKF-RS method by the 
relation 


Wens = U'^-W-UGE^''^^'=. (71) 

Threfore when the size of the ensemble is increased (by adding real or 
artihcial members) more information from the matrix W is captured by its 
projection onto the A^fc-dimensional space. Note that when —)■ n and 
U is orthonormal we have that Wens —t W. Consequently, the number 
of artihcial members will play an important role in the performance of the 
EnKF-RS implementation. 

4 Experimental Results 

This section tests the new EnKF implementations on a data assimilation 
problem using the quasi-geostrophic model presented in [38] . A comparison is 
done in two steps: hrst, the proposed implementations are compared against 
the well-known EnKF implementations presented in section and next the 
quality of the results for the EnKF-FS and EnKF-RS based on different 
values of N and K are assessed. 
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The oceans form a complex flow system influenced by the rotation of the 
Earth, the density stratihcation due to temperature and salinity, as well as 
other factors. The quasi-geostrophic (QG) model is a simple approximation 
of the real behavior of the ocean. It is dehned by the following partial 
differential equation: 

ojt + r ■ J {u, tjj) + (3 ■ — V ■ VV = —/r ■ + r ■ sin (”^2) 

in G G [0, Lx] X [0, Ly], where x and y represent the horizontal and vertical 
space components, u is the vorticity, -0 is the stream function, J (-0, u) is the 
Jacobian of two helds 


J (-0, u) = iIjx ■ Uy - tl^y ■ Ux , (73) 

and is the Laplacian operator. The coefficients 0, v, y and r are asso¬ 
ciated with the horizontal vorticity, the horizontal friction, the biharmonic 
horizontal friction, and the horizontal wind stress at the surface of the ocean, 
respectively. Moreover, the vorticity is related to the stream function by the 
elliptical equation: 


= uj. (74) 

The spatial domain for our experiments is G = [0, 1] x [0, 1]. The interior is 
covered by computational grids of different resolutions, and we denote by Di 
and D 2 the number of horizontal and vertical grid points, respectively. The 
different model resolutions are presented in the table 


Instance 

Di 

D 2 

n = Di- D 2 

<5^33x33 

31 

31 

961 

<5^65x65 

63 

63 

3,969 

QGi29x129 

127 

127 

16,129 


Table 1: Quasi-geostrophic instances for the computational tests in terms of 
the number of horizontal Di and vertical D 2 grid points. Di and D 2 do not 
consider boundary points. 

The numerical data assimilation experiments are characterized by the 
following settings: 
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Initial vorticities have the form 


ojq = sin(4 Xi Hj) ■ cos(2 Xi yj) + sin(2 Xi yj) + cos(4 Xi yj), 
ioT 1 < i < Di and 1 < j < D 2 . 

• The initial background error covariance is 

B = ■ Inxn, 

where the standard deviation a® is chosen to be 0.05 or 0.15 (times the 
true vorticity) in different experiments . 

• The observational errors are uncorrelated with variances [0.01]^. 

• The number of observed components from the vector state is given by 

m = p ■ n , (75) 

where p is the percentage of observed components. We consider two 
values for p, 0.7 and 0.9, corresponding to a sparser and a denser net¬ 
work, respectively. 

Other aspects of the numerical simulation are described below: 


• The EnKF methods are implemented in C-I--I- while the forward model 
(QG) [37] is implemented making use of FORTRAN. 

• The partial derivatives are discretized by central hnite differences. 


• The matrix and vector computations are efficiently carried out using 
the BLAS library |5]. 


• Matrix decompositions as well as eigenvalue computations are per¬ 
formed using the LAPACK library [1]. 


• The Arakawa method [2S] is utilized in order to compute the Jacobian 

(Pj). 


The time discretization of the model (72) uses of a fourth order Runge- 
Kutta method. The time step size 1.27 (units) which represent one 
hour in the ocean. The integration is performed for 1000 hours. 
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The NLOPT library [25] is utilized to numerically solve the optimiza¬ 
tion problems (39) and (43). 


The GSL-GNU scientihc library 
background and data errors. 


is utilized to generate the synthetic 


4.1 Comparison with current EnKF implementations 

We compare the EnKF-RS and EnKF-FS methods against current well- 
known EnKF implementations. The root mean square errors (RMSE) and 
GPU times for different methods making use of the QG^sxss, QGesxes and 
QGi 29 x 129 model instances are reported inTables|^[^ and|^ respectively. For 
a given number of snapshots and a reference trajectory 

where M is the number of assimilation times, the RMSE is dehned as follows 


RMSE = 




1 

M 


M 

E 

2=1 


X, 


^true I 


We vary the size of the ensemble N, the percentage of observed compo¬ 
nents p, and the initial background error cr®. As expected the traditional en¬ 
semble implementations EnKF, EnSRF, and EnTKF provide accurate anal¬ 
yses for different model resolutions and number of observed components. 
Moreover it can be seen in hgures and that the RMSE decreases as the 
simulation progresses. The traditional EnKF implementations provide the 
lowest elapsed time among the compared methods (i.e., EnKF and EnSRF 
implementations) which explains why they are attractive for use in real ap¬ 
plications. The RMSE values for the EnSRF and EnTKF implementations 
are identical since both hlters are deterministic and EnTKF is just an effi¬ 
cient implementation of the EnSRF. The inflation-free methods such as the 
EnKF-FN and EnK-DU implementations provide slightly better accurate re¬ 
sults than other methods. For example, for the largest instance < 5 ^ 129 x 129 , 
table 1^ shows the EnKF-DU to perform better than current implementa¬ 
tions. In most of the cases, strong duality holds: 


= VZs{C) and 


the slight differences between the optimal cost function values (39) and (43) 


are consistent with the numerical approximation errors in the solution of the 


optimization problems (40) and (45), respectively. 


Since the analysis state in the EnKF-DU formulation is obtained via the 


solution of the one-dimensional optimization problem (45), we expect this 
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method to be faster than the EnKF-FN implementation where the analysis 
requires the solution of the iV- dimensional optimization problem (39). This 
fact is also pointed out by Boquet in |6], Section 2.2], where the cost of com¬ 
puting the inverse (44) is assumed to be negligible in the dual formulation. 
However, this statement seems to be true for small model resolutions (i.e., 
Bocquet makes use of the Lorenz 96 model with 40 variables) and it holds for 
the smallest QG instance Nevertheless, for the QGgsxes case, the 

difference between the CPU times for the primal and dual implementations 
is almost negligible and even more, for the largest instance QG 129 X 129 , the 
EnKF-FN performs better than its dual approach for = 0.15. Although 
the cost function (43) depends only on (, every step in the optimization pro¬ 
cess requires the solution of the linear system (44) whose computational cost 
is not negligible in practice. Furthermore, when the initial background error 
is large the EnKF-DU computes many times the inverse of (44) and therefore 
its performance decreases considerably. 

Even if the traditional implementations perform very well in terms of 
RMSE and elapsed time, the most accurate results are obtained by the pro¬ 
posed new EnKF implementations. The results presented in hgures and [3^ 
show that the EnKF-RS method performs much better than traditional and 
inflation-free methods for the small instance QG^^xss, and are slightly better 
in the larger instances. The most accurate results among all the compared 
methods are the ones obtained by the EnKF-FS implementation. The results 
reported in hgures]^ andshow that, for all the instances and conhgurations, 
the RMSE obtained via the EnKF-FS outperforms the other methods by at 
least 60%. The CPU times for the proposed implementations are just slightly 
larger than those from the compared EnKF implementations. Since most of 
the computational time is spent in propagating the ensemble members, a 
modest increase in analysis time retains the potential of the new methods to 
perform well in practical applications. 
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RMSE 

CPU Time 

a 

B 

a 

B 

N 

P 

Method 

0.05 

0.15 

0.05 

0.15 



EnKF 

1.287 

3.858 

0.013 

0.009 



EnSRF 

1.289 

3.864 

0.013 

0.013 



EnTKF 

1.289 

3.864 

0.026 

0.028 


0.7 

EnKF-FN 

1.286 

3.855 

0.125 

0.116 



EnKF-DU 

1.285 

3.854 

0.045 

0.036 



EnKF-FS 

0.659 

1.974 

0.024 

0.036 

40 


EnKF-RS 

1.16 

3.49 

0.421 

0.252 


EnKF 

1.279 

3.836 

0.019 

0.021 



EnSRF 

1.282 

3.841 

0.01 

0.017 



EnTKF 

1.282 

3.841 

0.028 

0.017 


0.9 

EnKF-FN 

1.279 

3.831 

0.115 

0.115 



EnKF-DU 

1.276 

3.828 

0.072 

0.087 



EnKF-FS 

0.371 

1.116 

0.048 

0.028 



EnKF-RS 

1.07 

3.209 

0.395 

0.222 



EnKF 

1.268 

3.803 

0.065 

0.061 



EnSRF 

1.275 

3.82 

0.034 

0.052 



EnTKF 

1.275 

3.82 

0.143 

0.076 


0.7 

EnKF-FN 

1.264 

3.79 

0.608 

0.681 



EnKF-DU 

1.263 

3.79 

0.162 

0.181 



EnKF-FS 

0.646 

1.927 

0.119 

0.105 

80 


EnKF-RS 

1.069 

3.144 

1.303 

1.285 


EnKF 

1.252 

3.756 

0.074 

0.077 



EnSRF 

1.26 

3.773 

0.058 

0.051 



EnTKF 

1.26 

3.773 

0.168 

0.12 


0.9 

EnKF-FN 

1.249 

3.737 

0.894 

0.766 



EnKF-DU 

1.246 

3.737 

0.155 

0.181 



EnKF-FS 

0.358 

1.074 

0.192 

0.14 



EnKF-RS 

0.879 

2.588 

1.269 

1.26 


Table 2; RMSE and CPU-time (TIME) for the EnKF, EnSRF, EnTKF, 
EnKF-FN, EnKF-DU, EnKF-RS and EnKF-FS implementations applied to 
the QGssxsa instance (n = 961). 
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RMSE 

CPU Time 

a 

B 

a 

B 

N 

P 

Method 

0.05 

0.15 

0.05 

0.15 



EnKF 

1.67 

5.011 

0.072 

0.076 



EnSRF 

1.675 

5.025 

0.059 

0.047 



EnTKF 

1.675 

5.025 

0.056 

0.098 


0.7 

EnKF-FN 

1.671 

4.987 

0.185 

0.261 



EnKF-DU 

1.662 

4.987 

0.17 

0.245 



EnKF-FS 

0.893 

2.668 

0.108 

0.205 

40 


EnKF-RS 

1.61 

4.818 

1.869 

1.801 


EnKF 

1.667 

5.001 

0.095 

0.061 



EnSRF 

1.671 

5.015 

0.06 

0.035 



EnTKF 

1.671 

5.015 

0.067 

0.117 


0.9 

EnKF-FN 

1.667 

4.975 

0.262 

0.269 



EnKF-DU 

1.658 

4.975 

0.273 

0.244 



EnKF-FS 

0.52 

1.533 

0.209 

0.139 



EnKF-RS 

1.589 

4.761 

1.975 

2.011 



EnKF 

1.653 

4.958 

0.189 

0.189 



EnSRF 

1.661 

4.982 

0.174 

0.183 



EnTKF 

1.661 

4.982 

0.265 

0.377 


0.7 

EnKF-FN 

1.644 

4.915 

0.934 

0.919 



EnKF-DU 

1.638 

4.916 

0.567 

0.787 



EnKF-FS 

0.866 

2.592 

0.634 

0.495 

80 


EnKF-RS 

1.546 

4.612 

7.627 

7.542 


EnKF 

1.65 

4.949 

0.271 

0.302 



EnSRF 

1.657 

4.97 

0.227 

0.181 



EnTKF 

1.657 

4.97 

0.428 

0.429 


0.9 

EnKF-FN 

1.639 

4.888 

1.099 

1.01 



EnKF-DU 

1.63 

4.888 

0.85 

0.936 



EnKF-FS 

0.494 

1.486 

0.607 

0.502 



EnKF-RS 

1.507 

4.509 

7.845 

8.012 


Table 3: RMSE and CPU-time (TIME) for the EnKF, EnSRF, EnTKF, 
EnKF-FN, EnKF-DU, EnKF-RS and EnKF-FS implementations applied to 
the QGq 5 x 65 instance (n = 3969). 




RMSE 

CPU Time 

a 

B 

a 

B 

N 

P 

Method 

0.05 

0.15 

0.05 

0.15 



EnKF 

1.708 

5.125 

0.307 

0.186 



EnSRF 

1.712 

5.136 

0.174 

0.153 



EnTKF 

1.712 

5.136 

0.237 

0.223 


0.7 

EnKF-FN 

1.707 

5.1 

0.793 

0.589 



EnKF-DU 

1.7 

5.1 

0.788 

0.746 



EnKF-FS 

0.963 

2.847 

0.572 

0.611 

40 


EnKF-RS 

1.682 

5.042 

5.607 

5.832 


EnKF 

1.709 

5.127 

0.263 

0.277 



EnSRF 

1.712 

5.134 

0.166 

0.192 



EnTKF 

1.712 

5.134 

0.294 

0.296 


0.9 

EnKF-FN 

1.704 

5.089 

0.888 

0.849 



EnKF-DU 

1.696 

5.089 

1.047 

1.041 



EnKF-FS 

0.582 

1.661 

0.627 

0.745 



EnKF-RS 

1.676 

5.02 

6.368 

6.292 



EnKF 

1.702 

5.104 

0.721 

0.693 



EnSRF 

1.709 

5.125 

0.554 

0.576 



EnTKF 

1.709 

5.125 

0.838 

0.734 


0.7 

EnKF-FN 

1.69 

5.063 

3.085 

2.596 



EnKF-DU 

1.687 

5.063 

2.826 

2.931 



EnKF-FS 

0.941 

2.81 

2.191 

2.373 

80 


EnKF-RS 

1.654 

4.94 

26.384 

26.355 


EnKF 

1.7 

5.099 

0.834 

0.991 



EnSRF 

1.708 

5.123 

0.667 

0.729 



EnTKF 

1.708 

5.123 

0.924 

0.963 


0.9 

EnKF-FN 

1.683 

5.034 

4.07 

3.485 



EnKF-DU 

1.678 

5.034 

3.565 

3.501 



EnKF-FS 

0.558 

1.622 

2.745 

2.442 



EnKF-RS 

1.634 

4.908 

27.672 

26.208 


Table 4: RMSE and CPU-time for for the EnKF, EnSRF, EnTKF, EnKF- 
FN, EnKF-DU, EnKF-RS and EnKF-FS implementations applied to the 
<5^129x129 instance (n = 16129). 
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Figure 2: Plots of RMSE values for the EnKF, EnSRF, EnTKF, EnKF-FN, 
EnKF-DU, EnKF-FS and EnKF-RS for iV = 40 and p = 0.7. 
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Figure 3: Plots of RMSE values of the EnKF, EnSRF, EuTKF, EnKF-FN, 
EnKF-DU, EnKF-FS and EnKF-RS for iV = 80 and p = 0.7. 
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4.2 The impact of the number of ensemble members 
on performance 


The EnKF-RS and EnKF-FS implementations depend on the K samples 
taken from the distribution (57). We now study how the performance of the 
proposed methods varies for different values of N (number of real members) 
and K (number of artihcial members). For this we let K and N to be related 
by 


K = C-N, 


where C is a constant. Practical ensemble sizes range in 40 < iV < 80 
|22j . For each ensemble size we use several values of C between 0 and 10, 
e.g., = 40 and (P = 10 lead to K = 400. When C = 0 no artihcial 

members are added, but the error covariance matrix B is estimated. In 
the numerical experiments, the variances of the initial background error are 
set to = 0.15. The analyses RMSE values and the compute times for 
the proposed implementations using the QGasxss instance are reported in 
hgures and respectively. The analysis times for both implementations 
are small. Moreover, as expected, EnKF-RS is sensitive to changes in any 
of the parameters N and K. The RMSE is decreased when the values of 
those parameters are high as shown in hgures 4b and 4d The RMSE of the 


EnKF-FS analysis decreases only with increasing N as can be seen in hgures 
lialandlicl 

An important question is how well do the proposed implementations per¬ 
form with a small number of real members. Hopefully the inexpensive addi¬ 
tion of artihcial members can compensate for a small number of real ones. 
To this end we consider a small number of real members 10 < iV < 30 and 
a large number of artihcial members with 10 < C* < 60. The results for 


the EnKF-RS are shown in hgures ^ and for p equal to 0.7 and 0.9, re¬ 
spectively. The EnKF-RS implementation improves the estimated analysis 
state whenever N or K are increased. Moreover, the quality of the analyses 
obtained with small real ensembles (10 < N < 30) is comparable to those 
obtained with large real ensemble sizes (40 < N < 80). This justihes the 
addition of inexpensive artihcial members in order to increase the degrees 
of freedom of the ensemble. The EnKF-FS analysis improves only when the 
number of real members is increased. This can be seen in hgures 6a and 
6c For the smallest ensemble size {N = 10) the results obtained by the 
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EnKF-FS are better than those of any other implementation, inclnding the 
traditional EnKF implementations with the large ensemble size N = 80. 

Fignres ^ and M show that the EnKF-RS analyses for a small nnmber 
of observed components (70%) and a large nnmber of artihcial members 
[K) are eqnivalent to those obtained with a large nnmber of real members 
[N) and many observed components (90%). This is another indication of the 
positive impact obtained by increasing the nnmber of degrees of freedom with 
samples from the distribntion (57). This is compntationally less expensive 
than adding real members via rnnning the model. 



(c) EnKF-FS p = 0.9 (d) EnKF-RS p = 0.9 

Fignre 4: RMSE of the EnKF-FS and EnKF-RS implementations for different 
valnes of 0 < C < 10 and 40 < < 80 making nse of the QG 33 X 33 instance. 
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(c) EnKF-FS, p = 0.9 (d) EnKF-RS, p = 0.9 

Figure 5: Assimilation times of the EnKF-FS and EnKF-RS implementations 
for different values of 0 < C < 10 and 40 < < 80 making use of the 

QG33 x 33 instance 


34 

















(c) EnKF-FS, p = 0.9 (d) EnKF-RS, p = 0.9 


Figure 6 : RMSEs of the EnKF-FS and EnKF-RS implementations for dif¬ 
ferent values of 10 < (7 < 60 and 10 < iV < 30 making use of the QG 33 X 33 
instance. 
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(c) EnKF-FS, p = 0.9 (d) EnKF-RS, p = 0.9 

Figure 7: Assimilation times of the EnKF-FS and EnKF-RS implementations 
for different values of 10 < C* < 60 and 10 < < 30 making use of the 

QG'ss x 33 instance. 


5 Conclusions 

This paper develops two new implementations of the ensemble Kalman hlter 
(EnKF) based on shrinkage covariance estimation. The background error co- 
variance matrices used in analysis are obtained via the Rao-Blackwell Ledoit 
and Wolf estimator, which has been proved optimal in the estimation of 
high-dimensional covariance matrices from a small number of samples. This 
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covariance matrix and the background state (ensemble mean) serve as param¬ 
eters of the normal error distribution associated with the ensemble members. 
Samples from this distribution are taken in order to increase the number 
of ensemble members, and therefore, to decrease the sampling error in rep¬ 
resenting the background error distribution, and to increase the number of 
degrees of freedom in the assimilation. The two proposed implementations 
differ in the space where the assimilation process is performed: EnKF Full- 
Space (FnKF-FS) performs the analysis in the model space, while FnKF 
Reduce-Space (FnKF-RS) computes the analysis state in the space spanned 
by the ensemble members. Numerical experiments are carried out using 
a quasi-geostrophic model. They show that the two new implementations 
perform better than current FnKF implementations such as the traditional 
FnKF, square root hlters, and inflation-free FnKF methods. For all the 
scenarios and experimental settings, the FnKF-FS outperforms the other 
implementations by at least 60% in terms of accuracy (root mean square 
error). Moreover, for a small number of ensemble members (~ 10) and a 
moderate percentage of observed components from the vector state (~ 70%), 
the solutions obtained by the proposed methods are similar to those obtained 
by large ensemble sizes (~ 80) and large percentage of observed components 
(~ 90%). The computational time for analysis of the proposed implemen¬ 
tations is reasonably low. Since the total compute time is dominated by 
the multiple model runs considerable savings are expected from reducing the 
number of real ensemble members without deteriorating the quality of the 
results. 
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